classdef EddyDrawClass < handle
    % 用于Eddy实验过程出图
    %   此处显示详细说明
    
    properties
    % 研究区域
    extent
    % 研究区域的采样间隔
    interval
    filepath = 'F:\按天气候态\';
    end
    
    methods
        function obj = EddyDrawClass(varargin)
            %EDDYDRAWCLASS 构造此类的实例
            %   此处显示详细说明
            idx_extent = find(strcmp(varargin,'extent'));
            if ~isempty(idx_extent)
                obj.extent = varargin{idx_extent + 1};
            end

            idx_interval = find(strcmp(varargin,'interval'));
            if ~isempty(idx_interval)
                obj.interval = varargin{idx_interval + 1};
            end     
            idx_filepath = find(strcmp(varargin,'filepath'));
            if ~isempty(idx_filepath)
                obj.filepath = varargin{idx_filepath + 1};
            end  
        end
        
        function extentSet(obj, extent)
            obj.extent = extent;
        end
        
        function intervalSet(obj, interval)
            obj.interval = interval;
        end
        
        function drawEddyWithSSTA(obj, eddys, d)
            [lon_ssta, lat_ssta ,~ ,ssta, ~] = EddyArrayMatchSSTAClass.ssta_select_by_eddy(obj.filepath, eddys(1));
            item.lon = lon_ssta;
            item.lat = lat_ssta;
            item.data = ssta(:,:,d);
 
            %%%%%%%
            % res = 0.083;
            % lambda = 16;
            % item.data = filt2(item.data, res, lambda, 'hp');
            %%%%%%%

            item.extent = obj.extent;
            item.eddys = eddys;
            draw(item)
            function draw(item)
                lon = item.lon;
                lat = item.lat;
                data = item.data;
                extent = item.extent;
                
                [x,y]=meshgrid(lon,lat);
                x = x';
                y = y';
                % set(gcf,'position',[100 50 1200 800])
                % axes('position',[0.05 0.05 0.9 0.9])
                m_proj('miller','lon',[extent(2,1) extent(2,2)],'lat',[extent(1,1) extent(1,2)])
                m_pcolor(x,y,data);
                hold on;
                for i = 1:length(item.eddys)

                    s = m_scatter(item.eddys(i).longitude, item.eddys(i).latitude, 'filled','k');
                    s.SizeData = 5;
                    hold on;
                    if item.eddys(i).type  == -1
                        m_plot(item.eddys(i).contour_longitude, item.eddys(i).contour_latitude, '--k','LineWidth', 2);
                    else
                        m_plot(item.eddys(i).contour_longitude, item.eddys(i).contour_latitude, 'k','LineWidth', 2);
                    end
                    hold on;
                end
                % m_contour(x,y,data,'black','LineWidth',0.5,'showtext','on','levellist',[-2 -1 -0.5 0  0.5 1 2]);
                % hold on;
                xlabel('Longitude','FontSize',14);
                ylabel('Latitude','FontSize',14);
                m_gshhs_i('patch',[.6 .6 .6]);
                m_grid;
                colormap(nclCM(311,50));
                clim([-2 2]);
                colorbar;
            end

        end
        
        function drawEddyTraj(obj, traj)

            item.extent = obj.extent;
            item.traj = traj;
            draw(item);
            function draw(item)
                set(gcf,'position',[100 50 1200 800])
                axes('position',[0.05 0.05 0.9 0.9])
                m_proj('miller','lon',[item.extent(2,1) item.extent(2,2)],'lat',[item.extent(1,1) item.extent(1,2)])
                for i = 1:length(item.traj)
                    longitude = item.traj(i).longitude;
                    latitude = item.traj(i).latitude;
                    m_plot(longitude, latitude,'r');
                    hold on;
                end
                % m_contour(x,y,data,'black','LineWidth',0.5,'showtext','on','levellist',[-2 -1 -0.5 0  0.5 1 2]);
                % hold on;
                xlabel('Longitude','FontSize',14);
                ylabel('Latitude','FontSize',14);
                m_gshhs_i('patch',[.6 .6 .6]);
                m_grid;
            end

        end

        function DepthProfileDraw(obj, SSTA, depth, varargin)
            %% 画涡旋深度剖面
            [m,n,~] = size(SSTA);
            m_half = (m - 1)/2 + 1;
            n_half = (n - 1)/2 + 1;
            m_R = (m - 1)/20;
            n_R = (n - 1)/20;
            figure;
            set(gcf,'Units','normalized','Position',[0.2875,0.11,0.5175,0.815]);
            tiledlayout(1,2);
            nexttile
            X_pou = squeeze(SSTA(:,m_half,:));
            

            draw(X_pou, m_R, depth, 'X',varargin)
            nexttile
            Y_pou = squeeze(SSTA(n_half,:,:));


            draw(Y_pou, n_R, depth, 'Y',varargin)
            cb = colorbar;
            cb.Layout.Tile = 'east';
            title(cb, '\theta''/°C', 'FontSize', 9, 'FontWeight', 'bold', 'VerticalAlignment', 'bottom');
            function draw(pou, R, depth, str, varargin)
                if length(varargin)==1
                    varargin = varargin{1};
                end

                k = 1/9*ones(3);  % 用于降噪的核。个人理解类似于9点平均
                pou = conv2( pou , k , 'same' );
                
                x = -R:0.1:R;
                d = depth;
                [xx,dd] = meshgrid(x,d);
                pou = pou';
                % xx = xx';
                % dd = dd';
                % pcolor(xx,dd,pou);
                % hold on;
                % v = [-1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.1];
                [c,h] = contourf(xx,dd,pou,'black');
                % set(h,'LineColor','none');
                
                set(h,'LineWidth',0.1)%控制线宽
                set(h,'Fill','on')%填充
                % set(h,'levelstep',0.2)%线间距
                set(h,'showtext','on')%等高线标签
                xlabel(['\Delta ' str],'FontSize',14);
                ylabel('Depth','FontSize',14);

                % colormap(nclCM(310,200));%%
                colormap(nclCM(136,50));
                idx_clim = find(strcmp(varargin,'clim'));
                if ~isempty(idx_clim)
                    climExtent = varargin{idx_clim + 1};
                else
                    climExtent = [-2 2];
                end  

                clim(climExtent);
                xticks([-R:1:R]);
                yticks([0:100:1000]);
                xlim([-R R]);
                ylim([0 1000]);
                shading interp;
                axis ij;
            end
        end

        function  CrossSectionDraw(obj, SSTA, depth, depth_idx)
            %% 画涡旋深度截面
            [m,n,p] = size(SSTA);
            % m_half = (m - 1)/2 + 1;
            % n_half = (n - 1)/2 + 1;
            m_R = (m - 1)/20;
            n_R = (n - 1)/20;


            % Array = [0 50 80 100 150 200 300 500 700 1000];
            figure;
            set(gcf,'Position',[100,50,1500,900]);
            if length(depth_idx) == 1
                tiledlayout(1,length(depth_idx));
            else
                tiledlayout(2,length(depth_idx)/2);
            end
            for i = 1:length(depth_idx)
                nexttile
                surface = SSTA(:,:,depth_idx(i));
                draw(surface, depth(depth_idx(i)), m_R, n_R);
            
            end
            cb = colorbar;
            cb.Layout.Tile = 'east';
            title(cb, '\theta''/°C', 'FontSize', 9, 'FontWeight', 'bold', 'VerticalAlignment', 'bottom');
            
            function draw(surface, Depth, m_R, n_R)
               k = 1/25*ones(5);  % 用于降噪的核。个人理解类似于9点平均
                surface = conv2( surface , k , 'same' );
                x = -m_R:0.1:m_R;
                y = -n_R:0.1:n_R;
                [xx yy] = meshgrid(x,y);
                surface= surface';
                [c,h] = contourf(xx,yy,surface,'black','ShowText','on','TextList',[-2 -1.6 -1.2 -0.8 -0.4 -0.2 0 0.2 0.4 0.8 1.2 1.6  2],'levelstep',0.2,'LineWidth',0.1);
                hold on;
                txt = [num2str(Depth), 'm'];
                title(txt,'FontSize',14);
                ax = gca;
                ax.TitleHorizontalAlignment = 'left';
                xlabel('\Delta X','FontSize',14);
                ylabel('\Delta Y','FontSize',14);
                colormap(nclCM(382,50));
                xticks([-m_R:1:m_R]);
                yticks([-n_R:1:n_R]);
                clim([-2 2]);
                % shading interp;
                axis equal;
                axis xy;
            end
        end
        
        % 区域温度深度截面
        function ThermalSectionDraw(obj, area_SSTA, depth, varargin)
            % load('Mask.mat');
            [m,n] = size(area_SSTA);
            d = 36;
            pou_Array = zeros(m,n,d);
            for i = 1:m
                for j = 1:n
                    SSTA = area_SSTA(i,j).SSTA;
                    try
                    % Two_R_SSTA = SSTA(11:51,11:51,:);
                    One_R_SSTA = SSTA(21:41,21:41,:);
                    catch
                        pou_Array(i,j,:) = nan;
                        continue;
                    end
                    R = area_SSTA(i,j).R;
                    [~, SSTA_section] = EddyThermalTransortClass.thermalMeasure(One_R_SSTA, R, depth);
                    pou_Array(i,j,:) = squeeze(SSTA_section);
                end
            end


            % 画图重采样，默认为0.01间隔
            idx_depth = find(strcmp(varargin,'depth'));
            if ~isempty(idx_depth)
                depth = varargin{idx_depth + 1};
            else
                depth = 1;
            end  
            
            % figure;
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%
            surface = pou_Array(:,:,depth);

             
            % 画图重采样，默认为0.01间隔
            idx_reshapeInterval = find(strcmp(varargin,'reshapeInterval'));
            if ~isempty(idx_reshapeInterval)
                reshapeInterval = varargin{idx_reshapeInterval + 1};
            else
                reshapeInterval = 0.25;
            end   
            idx_boundary = find(strcmp(varargin,'boundary'));
            if ~isempty(idx_boundary)
                X = varargin{idx_boundary + 1};
                Y = varargin{idx_boundary + 1};
                mask = EddyThermalTransortClass.createMask(obj.extent,reshapeInterval,varargin);
                varargin{length(varargin) + 1} = 'mask';
                varargin{length(varargin) + 1} = mask;
            end   
            

            [~, ~, fi] = EddyThermalTransortClass.divand_interp(obj.extent, surface, reshapeInterval, varargin);

            lon = obj.extent(2,1) + 1/2 * reshapeInterval:reshapeInterval:obj.extent(2,2) - 1/2 * reshapeInterval;
            lat = obj.extent(1,1) + 1/2 * reshapeInterval:reshapeInterval:obj.extent(1,2) - 1/2 * reshapeInterval;

            item.lon = lon;
            item.lat = lat;
            item.data = fi;
            
            item.clim = [-2 2.10];
            item.extent = [min(item.lat) max(item.lat);min(item.lon) max(item.lon) ];
            draw2D(item, varargin);
            % set(gca,"FontName",'Times New Roman','FontSize',12);
            % cb = colorbar;
            % title(cb, '\theta''/°C', 'FontSize', 9, 'FontWeight', 'bold', 'VerticalAlignment', 'bottom');

            % saveas(f,'1.png');
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            function draw2D(item, varargin)

                %%%%%%% 南海九段线
                %中国省界线文件（含九段线）
                ChinaL=shaperead('H:\code\大论文数据处理\小论文画图\南海九段图\全国省级、地市级、县市级行政区划shp\全国省级、地市级、县市级行政区划shp\bou2_4l.shp');
                %中国面文件
                ChinaP=shaperead('H:\code\大论文数据处理\小论文画图\南海九段图\全国省级、地市级、县市级行政区划shp\全国省级、地市级、县市级行政区划shp\bou2_4p.shp');
                bou2_4lx=[ChinaL(:).X];
                bou2_4ly=[ChinaL(:).Y];
                bou2_4px=[ChinaP(:).X];
                bou2_4py=[ChinaP(:).Y];
                
                
                %%%%%%%

                if length(varargin)==1
                    varargin = varargin{1};
                end
                lon = item.lon;
                lat = item.lat;
                data = item.data;
                idx_picextent = find(strcmp(varargin,'picextent'));
                if ~isempty(idx_picextent)
                    extent = varargin{idx_picextent + 1};
                else
                    extent = item.extent;
                end
               
                data = smooth2a(data,5);
                [x,y]=meshgrid(lon,lat);
                x = x';
                y = y';
                % set(gcf,'position',[100 50 600 700])
                
                m_proj('miller','lon',[extent(2,1) extent(2,2)],'lat',[extent(1,1) extent(1,2)])
                % m_pcolor(x,y,data);
                [c,h] = m_contourf(x,y,data,'TextList',[-2 -1.6 -1.2 -0.8 -0.4 0  0.4 0.8 1.2 1.6 2],'ShowText','on');
                clabel(c,h,'FontSize',18,"FontName",'Times New Roman');
                % m_contour(x,y,data,'black','LineWidth',0.5,'showtext','on','levellist',[-2 -1 -0.5 0  0.5 1 2]);
                % hold on;
                set(h,'levelstep',0.2)%线间距
                % xlabel('Longitude','FontSize',14);
                % ylabel('Latitude','FontSize',14);
                % idx_title = find(strcmp(varargin,'title'));
                % if ~isempty(idx_title)
                %     titleName = varargin{idx_title + 1};
                %     title(titleName);
                % end   
                hold on;
                m_plot(bou2_4lx,bou2_4ly,'k','linewidth',1.5);%绘制中国省界
                clim(item.clim);

                m_gshhs_i('patch',[.6 .6 .6]);
                % m_grid('fontn','Times New Roman','fontsize',24);
             m_grid('fonts',24,'fontn','Times New Roman','ytick',[4 8 12 16 20 24], 'yticklabel',{"4° ","8° ","12° ","16° ","20° ","24°N"},'xtick',[108 111 114 117 120], 'xticklabel',{"108°","111°","114°","117°","120°E"});
                colormap(nclCM(231,50));
                % colormap(nclCM(149,200));
                % colorbar;
            %     saveas(f,saveName)
            %     close(f);
            end
        end

        function EHTDraw(obj, MHT, ZHT, varargin)
            
            MHT(find(isnan(MHT)))=0;
            ZHT(find(isnan(ZHT)))=0;
            
            MHT=medfilt2(MHT,[3 3]);
            ZHT=medfilt2(ZHT,[3 3]);
            
            
            lat = obj.extent(1,1):obj.interval:obj.extent(1,2);
            lon = obj.extent(2,1):obj.interval:obj.extent(2,2);

            % 画图重采样，默认为0.01间隔
            idx_reshapeInterval = find(strcmp(varargin,'reshapeInterval'));
            if ~isempty(idx_reshapeInterval)
                reshapeInterval = varargin{idx_reshapeInterval + 1};
            else
                reshapeInterval = 0.01;
            end   

            new_lat = obj.extent(1,1):reshapeInterval:obj.extent(1,2);
            new_lon = obj.extent(2,1):reshapeInterval:obj.extent(2,2);
            

            [x y] = meshgrid(lat(1:end-1),lon(1:end-1));
            [xx yy] = meshgrid(new_lat,new_lon);
            new_MHT = interp2(x,y,MHT,xx,yy);
            new_ZHT = interp2(x,y,ZHT,xx,yy);
            
            new_MHT=medfilt2(new_MHT,[3 3]);
            new_ZHT=medfilt2(new_ZHT,[3 3]);
            
            
            item.lon = new_lon;
            item.lat = new_lat;
            item.data = new_MHT/(10^6);
            item.clim = [-20 20];
            item.extent = [min(item.lon) max(item.lon) min(item.lat) max(item.lat)];
            f = figure(1)
            draw2D(item)
            title('经向运输');
            % saveas(f,'1.png');
            item.lon = new_lon;
            item.lat = new_lat;
            item.data = new_ZHT/(10^6);
            item.clim = [-20 20];
            item.extent = [min(item.lon) max(item.lon) min(item.lat) max(item.lat)];
            f = figure(2)
            draw2D(item)
            title('纬向运输');
            % saveas(f,'2.png');
            
            
            item.lon = lon(1:end-1);
            item.lat = lat(1:end-1);
            item.u = MHT/(10^6);
            item.v = ZHT/(10^6);
            item.clim = [-20 20];
            item.extent = [min(item.lon) max(item.lon) min(item.lat) max(item.lat)];
            f = figure(3)
            
            draw2Dquiver(item)
            title('温度运输场');
            % saveas(f,'3.png');
            
            [m,n] = size(MHT);
            MHT_integrate = zeros(1,n);
            ZHT_integrate = zeros(1,m);
            
            for i = 1:m
                for j = 1:n
                    if ~isnan(ZHT(i,j))
                        ZHT_integrate(i) = ZHT_integrate(i) + ZHT(i,j) * 0.25 * 111000;
                    end
                end
            end
            r        = 3;
            sigma    = 1;
            
            
            f = figure(4)
            x = lon(1:end-1);
            ZHT_integrate = Gaussianfilter(r, sigma, ZHT_integrate);
            plot(x,ZHT_integrate/ (10^12));
            ylim([-30 30]);
            title('沿经向积分的纬向输送');
            
            
            % saveas(f,'4.png');
            for i = 1:n
                for j = 1:m
                    if ~isnan(MHT(j,i))
                        MHT_integrate(i) = MHT_integrate(i) + MHT(j,i) * 0.25 * 111000;
                    end
                end
            end
            
            f = figure(5)
            x = lat(1:end-1);
            MHT_integrate = Gaussianfilter(r, sigma, MHT_integrate);
            plot(x,MHT_integrate/ (10^12));
            ylim([-30 30]);
            title('沿纬度积分的经向输送');
            % saveas(f,'5.png');
            
            
            function draw2D(item)
                draw2D_lon = item.lon;
                draw2d_lat = item.lat;
                data = item.data;
                extent = item.extent;
            
                [draw2D_x,draw2D_y]=meshgrid(draw2D_lon,draw2d_lat);
                draw2D_x = draw2D_x';
                draw2D_y = draw2D_y';
                set(gcf,'position',[100 50 1200 800])
                axes('position',[0.05 0.05 0.9 0.9])
                m_proj('miller','lon',[extent(1) extent(2)],'lat',[extent(3) extent(4)])
                m_pcolor(draw2D_x,draw2D_y,data);
            
                clim(item.clim);
                m_gshhs_i('patch',[.6 .6 .6]);
                m_grid;
                colormap(nclCM(151,200));
                % colormap(nclCM(149,200));
                colorbar;
            %     saveas(f,saveName)
            %     close(f);
            end
            
            
            
            function draw2Dquiver(item)
                draw2Dquiver_lon = item.lon;
                draw2Dquiver_lat = item.lat;
                u = item.u;
                v = item.v;
                extent = item.extent;
            
                % [x,y]=meshgrid(lon,lat);
                % x = x';
                % y = y';
                set(gcf,'position',[100 50 1200 800])
                axes('position',[0.05 0.05 0.9 0.9])
                m_proj('miller','lon',[extent(1) extent(2)],'lat',[extent(3) extent(4)])
                
                draw2Dquiver_new_lon = 105:0.01:122;
                draw2Dquiver_new_lat = 0:0.01:25;
                [draw2Dquiver_x, draw2Dquiver_y] = meshgrid(draw2Dquiver_lat,draw2Dquiver_lon);
                [draw2Dquiver_xx, draw2Dquiver_yy] = meshgrid(draw2Dquiver_new_lat,draw2Dquiver_new_lon);
            
                new_u = interp2(draw2Dquiver_x,draw2Dquiver_y,u,draw2Dquiver_xx,draw2Dquiver_yy);
                new_v = interp2(draw2Dquiver_x,draw2Dquiver_y,v,draw2Dquiver_xx,draw2Dquiver_yy);
            
            
                [draw2Dquiver_x, draw2Dquiver_y] = meshgrid(draw2Dquiver_lon,draw2Dquiver_lat);
                [draw2Dquiver_xx, draw2Dquiver_yy] = meshgrid(draw2Dquiver_new_lon,draw2Dquiver_new_lat);
                draw2Dquiver_x = draw2Dquiver_x';
                draw2Dquiver_y = draw2Dquiver_y';
                draw2Dquiver_xx = draw2Dquiver_xx';
                draw2Dquiver_yy = draw2Dquiver_yy';
                s = sqrt(new_u.^2 + new_v.^2);
                m_pcolor(draw2Dquiver_xx,draw2Dquiver_yy,s);
                hold on;
                m_quiver(draw2Dquiver_x,draw2Dquiver_y,u,v,3,'black');
                clim([0 20]);
                m_gshhs_i('patch',[.6 .6 .6]);
                m_grid;
                colormap(nclCM(127,200));
                % colormap(nclCM(149,200));
                colorbar;
            %     saveas(f,saveName)
            %     close(f);
            end
            
            
            % 功能：对一维信号的高斯滤波，头尾r/2的信号不进行滤波
            % r     :高斯模板的大小推荐奇数
            % sigma :标准差
            % y     :需要进行高斯滤波的序列
            function y_filted = Gaussianfilter(r, sigma, y)
            
                % 生成一维高斯滤波模板
                GaussTemp = ones(1,r*2-1);
                for i=1 : r*2-1
                    GaussTemp(i) = exp(-(i-r)^2/(2*sigma^2))/(sigma*sqrt(2*pi));
                end
                
                % 高斯滤波
                y_filted = y;
                for i = r : length(y)-r+1
                    y_filted(i) = y(i-r+1 : i+r-1)*GaussTemp';
                end
            end

            
        end
    end
end

